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A theory is developed for calculating vertical tunneling current between two sheets of bilayer graphene sep¬ 
arated by a thin, insulating layer of hexagonal boron nitride, neglecting many-body effects. Results are pre¬ 
sented using physical parameters that enable comparison of the theory with recently reported experimental 
results. Observed resonant tunneling and negative differential resistance in the current-voltage characteristics 
are explained in terms of the electrostatically-induced band gap, gate voltage modulation, density of states 
near the band edge, and resonances with the upper sub-band. These observations are compared to ones from 
similar heterostructures formed with monolayer graphene. 


In contrast to the well-known linear dispersion of 
monolayer graphene (MLG), charge carriers near the six 
corners of the Brillouin zone in an isolated graphene bi¬ 
layer are described by a quadratic energy dispersion. 1,2 
An even more intriguing distinction with MLG is that, 
under the influence of external fields, the band struc¬ 
ture of bilayer graphene (BLG) near the charge neutral¬ 
ity point becomes quartic, changing from semi-metallic 
to semiconducting as a small band gap is induced. 3-5 
This tunability of the band gap can be exploited by in¬ 
troducing gates, doping, and interactions with substrate 
materials in electronic devices based on BLG. 6-8 In this 
paper, we consider these effects and others in a 2D to 2D 
resonant tunneling device composed of two sheets of BLG 
separated by a thin, insulating layer of hexagonal boron 
nitride (h-BN). In a vertical configuration with an inter¬ 
layer bias, tunneling between two-dimensional electron 
gases is constrained by simultaneous energy and momen¬ 
tum conservation, leading to resonances in the current- 
voltage (/- V) characteristic and thus regions of negative 
differential resistance (NDR). 9 Such devices were orig¬ 
inally proposed for conventional 2D quantum wells, 10 
but the theory was recently treated for MLG, 11-18 and 
NDR was observed experimentally in high-quality de¬ 
vices shortly thereafter. 9,16,19 The theory discussed in the 
present work is particularly relevant to the recent obser¬ 
vations of Fallahazad et al. 9 

For a given interlayer voltage and for bilayers that are 
in crytallographic alignment, the electronic bands of the 
top and bottom bilayer of graphene will overlap for par¬ 
ticular sets of states with equal energy and crystal mo¬ 
mentum (Fig. 1). Away from the resonance voltage, only 
the states near the intersecting ring(s) can contribute to 
the tunneling current (Fig. lb). However, for one particu¬ 
lar voltage, the electrostatic potential between the bands 
will be zero, allowing all states between the two Fermi 
levels to tunnel simultaneously (Fig. lc). The shape and 
position of the resulting resonant peak in the I-V char¬ 
acteristic depends on the quantity and sign of charge car¬ 
riers in each bilayer, and therefore indirectly on external 
fields (gate voltages) and the electrostatic doping condi¬ 
tions. 

For example, in the absense of strong doping or sub- 



FIG. 1. (a) Device structure, with double black lines indicat¬ 
ing each graphene bilayer (BLG) and the group of orange lines 
representing 4 to 6 layers of h-BN. (b) Alignment of electronic 
bands at an off-resonant interlayer bias voltage; blue (solid) 
curves for one bilayer and red (dashed) for the other; (c) at 
the voltage which yields the primary tunneling resonance; (d) 
at a higher voltage which aligns the lower bands of one bilayer 
with the higher sub-bands of the other. The largest contri¬ 
bution to tunneling current occurs near the states where the 
two bands intersect. Bands represent energy as a function of 
in-plane crystal momentum near one of the six corners of the 
Brillouin zone. 


strate interactions, resonance can be observed for both 
both positive and negative bias voltages as the back-gate 
voltage (Vbg) is swept from one sign to the other (Fig. 2). 
Recently, Fallahazad et al. have observed resonances 
with precisely this behavior in devices fabricated with 
exfoliated BLG/h-BN/BLG on a h-BN/Si02 substrate. 9 
The width and amplitude of each resonant peak relative 
to the background (non-resonant) current are determined 
by the degree of coherence between tunneling wavefunc- 
tions, as is discussed in detail in Ref. 17. 
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FIG. 2. Calculated tunnel current density versus interlayer 
bias for undoped graphene bilayers for a range of gate volt¬ 
ages. Upper inset: a similar computation for a larger volt¬ 
age range highlighting secondary resonances from the higher 
sub-bands. Lower inset: a closer view of the Fbg = 0 case, 
varying coherence length, a disorder parameter in the theory, 
from 50nm (solid, light) to 10nm (solid, dark). 


In addition to the primary resonance, the higher sub¬ 
band of one bilayer can also come into alignment with 
the lower sub-band of the second bilayer causing a simi¬ 
lar spike in the tunneling current. Secondary resonances 
as well as an increase in background current from the up¬ 
per bands entering the tunneling energy window can be 
observed at larger voltages as shown in the upper inset of 
Fig. 2. Interactions with the upper sub-bands are distinct 
from monolayer devices, and may provide opportunities 
for multi-state logic devices. 

At a smaller voltage scale, and especially at lower tem¬ 
peratures, it is possible to observe additional features 
due to the tunable band gap in BLG. The presence of 
a transverse electric field across a graphene bilayer in¬ 
duces a potential difference between the two individual 
layers of graphene. This broken layer symmetry causes 
a small band gap to open up around the charge neu¬ 
trality point which increases with the magnitude of the 
potential difference across the bilayer. In the tunnel¬ 
ing device modeled here, the interlayer and gate volt¬ 
ages modulate the separate potential difference across 
each individual bilayer in a coupled system . 8 As a re¬ 
sult, the band gaps in both bilayers vary with voltage 
(typically at different rates), which affects the overall 
tunneling current. For non-zero band gap, the precise 
form of the energy dispersion is quartic near the gap, as 
in Figs. 1 and 3(a). Moreover, the location of the band 



ps(E) pd{E ) 


FIG. 3. (a) Electronic bands in the source (left) and drain 
(right) electrodes with the tunneling barrier (band gap of 
boron nitride) in between at a small positive bias. Dashed 
lines indicate the Fermi levels in each bilayer, /j, = —eld; not 
to scale. Density of states corresponding to (b) the source 
electrode and (c) the drain electrode in the same bias condi¬ 
tion as panel (a); energy axes in units of eV. The alignment 
of the divergences in the density of states near the valence 
band edge of each bilayer produces a large overlap of states 
and thus a spike in tunneling current. 


gap is a ring of states concentric with the AT-point. This 
arrangement of states causes divergences in the DOS at 
the conduction and valence band edges, which can yield 
additional spikes in the tunneling current for certain elec¬ 
trostatic arrangements. Whereas the primary feature in 
the tunneling current occurs when the electrostatic po¬ 
tentials in the source <ps and drain (f>u electrodes are 
aligned, A</> = (/)$ — 4>d = 0 , other features due to over¬ 
lap of the large DOS near the band edges can occur when 
one of the four conditions A cj) ± E 9y s /2 ± E g ^/2 = 0 is 
satisfied (where E g 7 is the band gap in each bilayer), as 
in Figs. 3 and 4. These features in the I-V characteristic 
are distinct from those caused directly by momentum- 
conserving effects with complete band alignment (as in 
Fig. 1) and are less sensitive to the relaxation of momen¬ 
tum conservation (decoherence), but may be observed in 
tandem with the latter. In MLG there are no equivalent 
band edges, and thus these additional sharp features are 
absent in monolayer vertical tunneling devices. 

We use a tight-binding model for the dispersion of BLG 
with nearest-neighbor hopping energy 70 « 3.1eV, inter¬ 
layer hopping energy 71 ss 0.4 eV, and interlayer poten¬ 
tial asymmetry U . 20 Higher order considerations such as 
the trigonal warping of the bands (azimuthal asymme¬ 
try) were found to have a negligible impact on the tun¬ 
neling and thus were excluded. The occupation of levels 
and band gap in each electrode varies with the set of ap¬ 
plied voltages, and thus the electrostatic potentials are 
required to calculate the tunneling current. These poten¬ 
tials are determined by first solving a matrix equation 
<7i = CijVj , treating each monolayer of graphene sepa¬ 
rately, to obtain the transverse fields across each bilayer. 
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FIG. 4. Low-voltage tunneling current for a device with a 
small amount of doping on the top (drain) bilayer at 10 K 
showing a number of small features due to the alignment of 
various band-edges, as explained in Fig. 3. I-V curves are 
shifted vertically for clarity. Numbered insets show the band 
alignment for each of the four labeled points along the Vbg = 
40 V curve. Arrows indicate electron current that produces 
the sharp feature in each case. 


We then use those fields to solve a second matrix equa¬ 
tion treating each bilayer with the local DOS for each 
layer within the bilayers. This method can accomodate 
both top and bottom gates, though we chose to focus on 
matching with devices with only one gate in the present 
work. Net charges are calculated using full Fermi inte¬ 
grals = e(n i -p J )-e(u 0) i-po,i), rii = e f d E p{E)f{E) 
to account for quantum capacitance and thermal occupa¬ 
tion, with environmental doping densities nop. We calcu¬ 
late the tunneling current by summing over the transition 
rates between all states in the source and drain bilayers, 

i = gsgv^Y, \ M <*e\ 2 \fs{E a ) - f D (E p )]S(E a - e 0 ), 

a,/3 

(i) 

with spin and valley degeneracies g s , g v and state labels 
a and (3 in the source and drain bilayers. 11 The over¬ 
lap integrals between states in the source and drain are 


contained in the matrix element 


M a p 


h 2 _ 

2m 


d S 




( 2 ) 


which is evaluated in a similar way as for MLG in Refs. 
11 and 17. We calculate the surface integral in Eq. 2 
over a region defined by the length scale of wavefunc- 
tion coherence in the device, a parameter we call the 
characteristic coherence length, L. This is a disorder 
parameter which defines the degree of momentum con¬ 
servation and thus controls the width and amplitude of 
resonant features in the I—V characteristic. The momen¬ 
tum (wavevector) conservation, chiral (angular) terms, 
and crystallographic misorientation are encapsulated in 
the matrix elements, while energy conservation is con¬ 
tained in the ^-function that appears in Eq. I. 11,1 ' 

In contrast to the theory for MLG devices, for BLG 
this ^-function must be evaluated using the quartic dis¬ 
persion relation in order to capture band-gap and higher 
sub-band effects. Converting the sums over states in Eq. 
1 to integrals over fc, we can evaluate the (5-function by 
changing variables from E to k, 

6[E(k s ) - E(k D ) + A0] ^ 


for all combinations of bands between each bilayer, where 
/ is equal to the original argument of the (5-function, ki 
are the zeros of this argument, and /' is the derivative 
with respect to k. This procedure allows us to remove 
one fc-integration and proceed to calculating the current. 
A small amount of broadening is introduced to handle the 
singularities that arise near the band edges (an imaginary 
term it is added to the \f'\ denominator, with epsilon 
typically equal to 10 ~ 2 Tivf)- 

Comparing our theory with the experimental results of 
Fallahazad et al. 9 , we find for the undoped device at room 
temperature (Fig. 2) very good agreement both in terms 
of the peak shapes and the gate-voltage dependence. For 
the low-temperature results of Fig. 4, small peaks as¬ 
sociated with DOS features become prominent, super¬ 
imposed on a broad momentum-conserving background 
current. We believe the situation found in experiment 
at low temperature is the same, showing a similar sharp 
peak superimposed on a smooth background current. 
The interpretation offered in Ref. 9 associates the sharp 
peak itself with a momentum-conserving resonant effect, 
but no origin for the broad background is provided. Al¬ 
ternatively, in our interpretation, both features can be 
well understood. The data for the undoped device at 
room temperature (Fig. 2) can be similarly understood 
within the same framework. Sharp DOS features are not 
seen for the latter, either in theory or experiment, since 
the higher temperature leads to a reduction in the am¬ 
plitude of the sharp peaks (at elevated temperature the 
tunnel current includes contributions from nearby states 
that are thermally occupied, leading to a reduction in 
strength of the sharp peaks). This distinction between 
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DOS versus momentum-conserving effects, as provided 
by our theory, provides an expanded interpretation of 
the experimental results. 9 

While resonant tunneling in MLG heterostructures is 
novel and intriguing, the additional sub-bands in BLG 
as well as its unusual behavior in the presence of trans¬ 
verse fields provides many additional channels for inter¬ 
esting tunneling phenomena. Although the results pre¬ 
sented here were calculated with zero angular misorien- 
tation (perfect crystallographic alignment) between the 
two bilayers of graphene, the theory readily computes 
current for non-zero misorientation, as discussed and ob¬ 
served in prior work for MLG. 11,17,19 Concerning possible 
misorientation within the graphene bilayers themselves, 
this is known experimentally not to occur for the devices 
of Fallahazad et al. 9 An additional source of misorienta¬ 
tion in the device would be that between the graphene 
bilayers and the h-BN layers of the tunnel barrier. We 
have not investigated this effect in detail, although refer¬ 
ring to prior work for twisted BLG, 21,22 it appears that 
such an effect would give rise to a reduced transmission 
current through the entire heterostructure. Indeed, for 
the case of tunneling between MLG layers separated by 
h-BN, computed tunnel currents agree in detail with ex¬ 
periment, except that the theory is a factor of 10 3 to 
10 4 too large. 17 We find a similar discrepancy in abso¬ 
lute magnitude of the current for the present situation 
of BLG/h-BN devices, and we consider it likely that the 
reduced conductance of the BLG/h-BN interface is the 
source of this discrepancy. 

For BLG devices, we find that DOS effects are largely 
unaffected by small amounts of angular misorientation 
between the bilayers, whereas momentum-conserving res¬ 
onant peaks are shifted due to the change in conditions 
required for band intersection, as in monolayer devices. 
We note that the electronic properties of the BLG can 
be expected to be influenced by the neighboring h-BN, 
in analogy with the MLG case. 23 Such effects are typi¬ 
cally on the 1 to lOmeV scale; they will be important 
for a very detailed comparison between experiment and 
theory, but in terms of the overall distinction that we 
make here between DOS and momentum-conserving ef¬ 
fects these effects can be neglected. Similarly, we neglect 
many-body modifications to the BLG band structure (in¬ 
cluding many-body effects between the two graphene bi¬ 
layers, since they are separated in the experiments 9 by 4 
to 6 monolayers of h-BN). The effects of external in-plane 
magnetic fields have been explored for similar monolayer 
and monolayer/bilayer devices, 18,19,24 but are not consid¬ 
ered here for brevity. Finally, inelastic effects may play a 
role in some devices, particularly at room temperature, 
however, we have focused here on elastic interactions, 
which play a large role in the relaxation of momentum 
conservation and subsequently the strength of resonant 
behavior compared to background current. 
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